//    Copyright (C) Kherty.  All rights reserved.
using System;

namespace OpenLS.Spreadsheet.StandardAddIn.Statistics
{
    internal class HyperGeometricDistribution : IntegerDistribution
    {
        private int _numberOfSuccesses;
        private int _populationSize;
        private int _sampleSize;

        public HyperGeometricDistribution(int populationSize, int numberOfSuccesses, int sampleSize)
        {
            if (numberOfSuccesses > populationSize)
                throw new ArgumentOutOfRangeException("Number of success must be less or equal to population size");
            if (sampleSize > populationSize)
                throw new ArgumentOutOfRangeException("Sample must be less than or equal to population size");
            PopulationSize = populationSize;
            SampleSize = sampleSize;
            NumberOfSuccesses = numberOfSuccesses;
        }

        public int NumberOfSuccesses
        {
            get { return _numberOfSuccesses; }
            set
            {
                if (value < 0)
                {
                    throw new ArgumentOutOfRangeException("value");
                }
                _numberOfSuccesses = value;
            }
        }

        public int PopulationSize
        {
            get { return _populationSize; }
            set
            {
                if (value <= 0)
                    throw new ArgumentOutOfRangeException("value");
                _populationSize = value;
            }
        }

        public int SampleSize
        {
            get { return _sampleSize; }
            set
            {
                if (value < 0)
                    throw new ArgumentOutOfRangeException("value");
                _sampleSize = value;
            }
        }

        public override double CumulativeProbability(int x)
        {
            int n = PopulationSize;
            int m = NumberOfSuccesses;
            int k = SampleSize;

            int[] domain = getDomain(n, m, k);
            if (x < domain[0])
                return 0;
            else if (x >= domain[1])
                return 1.0;
            return innerCumulativeProbability(domain[0], x, 1, n, m, k);
        }

        private static int[] getDomain(int n, int m, int k)
        {
            return new int[]
                {
                    GetLowerDomain(n, m, k),
                    GetUpperDomain(m, k)
                };
        }

        protected override int DomainLowerBound(double p)
        {
            return GetLowerDomain(PopulationSize, NumberOfSuccesses,
                                  SampleSize);
        }

        protected override int DomainUpperBound(double p)
        {
            return GetUpperDomain(SampleSize, NumberOfSuccesses);
        }

        private static int GetLowerDomain(int n, int m, int k)
        {
            return Math.Max(0, m - (n - k));
        }

        private static int GetUpperDomain(int m, int k)
        {
            return Math.Min(k, m);
        }

        public override double Probability(int x)
        {
            double result;

            int n = PopulationSize;
            int m = NumberOfSuccesses;
            int k = SampleSize;

            int[] domain = getDomain(n, m, k);
            if (x < domain[0] || x > domain[1])
            {
                result = 0.0;
            }
            else
            {
                result = probability(n, m, k, x);
            }

            return result;
        }

        private static double probability(int n, int m, int k, int x)
        {
            return Math.Exp(Statistician.BinomialCoefficientLog(m, x) +
                            Statistician.BinomialCoefficientLog(n - m, k - x) -
                            Statistician.BinomialCoefficientLog(n, k));
        }

        public double UpperCumulativeProbability(int x)
        {
            double result;

            int n = PopulationSize;
            int m = NumberOfSuccesses;
            int k = SampleSize;

            int[] domain = getDomain(n, m, k);
            if (x < domain[0])
            {
                result = 1.0;
            }
            else if (x > domain[1])
            {
                result = 0.0;
            }
            else
            {
                result = innerCumulativeProbability(domain[1], x, -1, n, m, k);
            }

            return result;
        }

        private double innerCumulativeProbability(
            int x0, int x1, int dx, int n, int m, int k)
        {
            double result = probability(n, m, k, x0);
            while (x0 != x1)
            {
                x0 += dx;
                result += probability(n, m, k, x0);
            }
            return result;
        }
    }
}